Evolution of phenotypic plasticity leads to tumor heterogeneity with implications for therapy

Cancer is a significant global health issue, with treatment challenges arising from intratumor heterogeneity. This heterogeneity stems mainly from somatic evolution, causing genetic diversity within the tumor, and phenotypic plasticity of tumor cells leading to reversible phenotypic changes. However, the interplay of both factors has not been rigorously investigated. Here, we examine the complex relationship between somatic evolution and phenotypic plasticity, explicitly focusing on the interplay between cell migration and proliferation. This type of phenotypic plasticity is essential in glioblastoma, the most aggressive form of brain tumor. We propose that somatic evolution alters the regulation of phenotypic plasticity in tumor cells, specifically the reaction to changes in the microenvironment. We study this hypothesis using a novel, spatially explicit model that tracks individual cells’ phenotypic and genetic states. We assume cells change between migratory and proliferative states controlled by inherited and mutation-driven genotypes and the cells’ microenvironment. We observe that cells at the tumor edge evolve to favor migration over proliferation and vice versa in the tumor bulk. Notably, different genetic configurations can result in this pattern of phenotypic heterogeneity. We analytically predict the outcome of the evolutionary process, showing that it depends on the tumor microenvironment. Synthetic tumors display varying levels of genetic and phenotypic heterogeneity, which we show are predictors of tumor recurrence time after treatment. Interestingly, higher phenotypic heterogeneity predicts poor treatment outcomes, unlike genetic heterogeneity. Our research offers a novel explanation for heterogeneous patterns of tumor recurrence in glioblastoma patients.


Overview
The model dynamics is defined via two operators: a deterministic transport operator T , and a stochastic rearrangement operator R. The latter comprises births and deaths, phenotypic switches, and changes in migration velocity.The transport operator propagates cells in velocity channels towards their respective direction of movement.
The composition of the two operators T • R is applied independently at every node r and at each time step k to create the configuration at the next time step k + 1.
The stochastic rearrangement operator is defined as a composition of four processes: cell death D, phenotypic switching S, proliferation B, and orientation O, The precise updated rules are defined below.Note that we assume that each cell dies, switches phenotype, proliferates, and switches direction of movement independently.

Cell death
We assume that each tumor cell dies independently from other cells with probability δ in each time step.Therefore, the number of dying cells with property p follows a binomial distribution, and we can write the new configuration after cell deaths as where η D is a configuration that contains each dying cell.The sum of channel occupations a + b, where a, b ∈ M P , is the multiset c ∈ M P defined by its multiplicity Similarly, the difference of two channel occupations a − b, where a, b ∈ M P , as c ∈ M P is defined as via the multiplicity as The definitions of sum and difference of channel occupations extend to configurations of nodes by point-wise application, i. e., c = a ± b, with a, b, c ∈ E, is defined as η D , the node configuration indicating dying cells, is formally defined by the probability where M ηi (p) indicates the number of cells with property p in channel i of node configuration η.

Phenotypic switch
After the death operator has removed the dying cells, cells switch their phenotype according to the phenotypic switch probability r κ (ρ N ) ∈ [0, 1] depending on the individual switch parameter κ and the average local density ρ N (r, k) := To express this process mathematically, we first define the multiset of all cells on a node Based on this, we construct the multisets of cells with the migratory η m and proliferative η p phenotypes after the phenotype switch, respectively

Proliferation
Cells in the proliferative phenotype produce offspring with probability The offspring is initially in the proliferative phenotype and contained in the offspring multiset η B .Each daughter cell's phenotypic switch parameter κ ′ is drawn from a Gaussian probability distribution f (κ ′ |κ) centered around their mother cell's phenotypic switch parameter κ κ ′ ∼ N (κ, ∆κ).(S12) The multiset of cells in the proliferative phenotype after the birth step is then The offspring multiset η B is informally defined as a stochastic process in the following way

Reorientation
In the final substep of the stochastic rearrangement, cells with the migratory phenotype choose their new direction of movement independently and with uniform probability T i = 1/b.Consequently, the probability for the post-orientation configuration η R is a multinomial distribution defined by the transition probability Here δ(x, y) := 1 if x = y, 0 else , is the Kronecker delta, which ensures mass conservation during the orientation step and must not to be confused with the death probability δ.The rest channel occupation is given by the multiset of cells with the proliferative phenotype after the birth step and unchanged by the orientation step, i. e., η R b = η ′′′ p .Therefore the second product in ?? is only multiplying the terms corresponding to the velocity channels (the index runs from i = 0 to i = b − 1).

Transport
After stochastic rearrangement, migratory cells move in their direction of movement, while proliferative cells stay in the rest channel of the same node.In other words

Derivation of per-capita growth rate
Here, we derive the average per-capita growth rate of a cell with switch parameter κ in a microenvironment with mean density ρ N .The probability for death is equal to δ, i. e., P (∆n = −1) = δ.The probability for proliferation is equal to the probability of survival times the probability of being in the proliferative phenotype times the probability of proliferation, i. e., P (∆n = +1) = (1 − δ)α(1 − ρ)r κ (ρ N ).Therefore, the average change is 3 Scaling of LGCA dimensions

Maximum tumor size and cell density
Here, we estimate the maximum tumor size of our synthetic tumors in physical units.
To do so, we treat our one-dimensional simulations as projections along a spherically growing tumor.We assume a typical glioblastoma radius of cancer cells of r cell ≈ 7 µm [1].This allows us to estimate the volume of one LGCA voxel via its capacity K, assuming that the voxel is fully occupied with cells at the carrying capacity.Therefore V node = KV cell = 100 • 4 3 πr 3 cell ≈ 1.4•10 5 µm 3 .Thus, we estimate the side length of one voxel as L node ≈ 50 µm.If the lattice is almost fully filled at the end of our simulations and we start with cells in the middle of the lattice, the maximum final tumor diameter is approximately L • L node = 1001 • 50 µm ≈ 5 cm.
These parameters correspond to a tumor cell density of K/V node ≈ 10 6 mm −3 , which is the same order of magnitude as in glioma patients [2].

Cell cycle and recurrence times
To estimate the time step length in real units, we consider the cell cycle time of glioblastoma cells in vitro.For these cells, effective cell cycle times (neglecting births of non-viable daughter cells) T b = 38±4 h were reported [3].This value should correspond to that of a single cell in our model, i. e., one that is not inhibited in its growth by other cells.The cell cycle time and growth rate α are connected using the time step , k) ∈ N 0 Cell number at node r and time step k ϵ > 0 Lattice spacing τ > 0 Time step length ρ N (r, k) > 0 Average density around node r at time step k δ ∈ [0, 1] Death probability α(1 − n(r, k)/K) Proliferation probability of resting cells r κ ∈ [0, 1] Probability of switching to the proliferative phenotype 1 − r κ ∈ [0, 1] Probability of switching to the migratory phenotype κ ∈ R Individual phenotypic switch regulation θ ∈ [0, 1] Critical cell density where phenotypes are equally probable ∆κ ∈ R + Standard deviation of κ-values of offspring versus mother cells T b > 0 Cell cycle time of glioblastoma cells in vitro length τ via α(τ ) = τ T b ⇔ τ = αT b .(S18) For α = 1 we have τ = T b ≈ 38 h ≈ 1.6 d.